In [2]:
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os
In [3]:
climate = pd.read_csv(r'C:\Users\kacie\Box\Personal\CareerFoundry\Machine Learning\Climate Data\Scaled Climate Data.csv')
In [4]:
climate.head()
Out[4]:
BASEL_cloud_cover BASEL_wind_speed BASEL_humidity BASEL_pressure BASEL_global_radiation BASEL_precipitation BASEL_snow_depth BASEL_sunshine BASEL_temp_mean BASEL_temp_min ... VALENTIA_pressure VALENTIA_global_radiation VALENTIA_precipitation VALENTIA_snow_depth VALENTIA_sunshine VALENTIA_temp_mean VALENTIA_temp_min VALENTIA_temp_max Date Month
0 0.660514 -0.02793 0.826097 -0.001949 -1.101066 -0.265148 -0.179228 -0.902918 -0.528623 -0.845652 ... -1.299744 -0.806427 -0.088407 -0.024706 0.372147 -0.668215 -0.519743 -0.752237 19600101 1
1 0.244897 -0.02793 0.735760 -0.001949 -1.058108 1.658760 -0.179228 -0.810126 -0.582946 -0.462450 ... -1.262455 -1.042055 0.503361 -0.024706 -0.829285 -0.548046 -0.629054 -0.407141 19600102 1
2 1.076130 -0.02793 1.277781 -0.001949 -1.251420 0.155707 -0.179228 -1.065304 -0.257010 -0.186545 ... -0.432779 -1.136306 -0.396127 -0.024706 -1.009500 -0.067372 0.054135 -0.177078 19600103 1
3 -1.001953 -0.02793 1.458455 -0.001949 -0.821838 -0.445514 -0.179228 -0.114186 -0.555784 -0.385810 ... 0.387574 -1.183432 0.669056 -0.024706 -1.039536 -0.998679 -0.164486 -0.838511 19600104 1
4 0.244897 -0.02793 1.729466 -0.001949 -0.746661 -0.164944 -0.179228 0.187388 -1.003946 -1.075573 ... 1.729970 -0.794645 -0.490810 -0.024706 0.672505 -1.509396 -1.339569 -1.471186 19600105 1

5 rows × 170 columns

In [5]:
# view column names to find temperature columns
print(climate.columns.to_list())

['BASEL_cloud_cover', 'BASEL_wind_speed', 'BASEL_humidity', 'BASEL_pressure', 'BASEL_global_radiation', 'BASEL_precipitation', 'BASEL_snow_depth', 'BASEL_sunshine', 'BASEL_temp_mean', 'BASEL_temp_min', 'BASEL_temp_max', 'BELGRADE_cloud_cover', 'BELGRADE_humidity', 'BELGRADE_pressure', 'BELGRADE_global_radiation', 'BELGRADE_precipitation', 'BELGRADE_sunshine', 'BELGRADE_temp_mean', 'BELGRADE_temp_min', 'BELGRADE_temp_max', 'BUDAPEST_cloud_cover', 'BUDAPEST_humidity', 'BUDAPEST_pressure', 'BUDAPEST_global_radiation', 'BUDAPEST_precipitation', 'BUDAPEST_sunshine', 'BUDAPEST_temp_mean', 'BUDAPEST_temp_min', 'BUDAPEST_temp_max', 'DEBILT_cloud_cover', 'DEBILT_wind_speed', 'DEBILT_humidity', 'DEBILT_pressure', 'DEBILT_global_radiation', 'DEBILT_precipitation', 'DEBILT_sunshine', 'DEBILT_temp_mean', 'DEBILT_temp_min', 'DEBILT_temp_max', 'DUSSELDORF_cloud_cover', 'DUSSELDORF_wind_speed', 'DUSSELDORF_humidity', 'DUSSELDORF_pressure', 'DUSSELDORF_global_radiation', 'DUSSELDORF_precipitation', 'DUSSELDORF_snow_depth', 'DUSSELDORF_sunshine', 'DUSSELDORF_temp_mean', 'DUSSELDORF_temp_min', 'DUSSELDORF_temp_max', 'GDANSK_cloud_cover', 'GDANSK_humidity', 'GDANSK_precipitation', 'GDANSK_snow_depth', 'GDANSK_temp_mean', 'GDANSK_temp_min', 'GDANSK_temp_max', 'HEATHROW_cloud_cover', 'HEATHROW_humidity', 'HEATHROW_pressure', 'HEATHROW_global_radiation', 'HEATHROW_precipitation', 'HEATHROW_snow_depth', 'HEATHROW_sunshine', 'HEATHROW_temp_mean', 'HEATHROW_temp_min', 'HEATHROW_temp_max', 'KASSEL_wind_speed', 'KASSEL_humidity', 'KASSEL_pressure', 'KASSEL_global_radiation', 'KASSEL_precipitation', 'KASSEL_sunshine', 'KASSEL_temp_mean', 'KASSEL_temp_min', 'KASSEL_temp_max', 'LJUBLJANA_cloud_cover', 'LJUBLJANA_wind_speed', 'LJUBLJANA_humidity', 'LJUBLJANA_pressure', 'LJUBLJANA_global_radiation', 'LJUBLJANA_precipitation', 'LJUBLJANA_sunshine', 'LJUBLJANA_temp_mean', 'LJUBLJANA_temp_min', 'LJUBLJANA_temp_max', 'MAASTRICHT_cloud_cover', 'MAASTRICHT_wind_speed', 'MAASTRICHT_humidity', 'MAASTRICHT_pressure', 'MAASTRICHT_global_radiation', 'MAASTRICHT_precipitation', 'MAASTRICHT_sunshine', 'MAASTRICHT_temp_mean', 'MAASTRICHT_temp_min', 'MAASTRICHT_temp_max', 'MADRID_cloud_cover', 'MADRID_wind_speed', 'MADRID_humidity', 'MADRID_pressure', 'MADRID_global_radiation', 'MADRID_precipitation', 'MADRID_sunshine', 'MADRID_temp_mean', 'MADRID_temp_min', 'MADRID_temp_max', 'MUNCHENB_cloud_cover', 'MUNCHENB_humidity', 'MUNCHENB_global_radiation', 'MUNCHENB_precipitation', 'MUNCHENB_snow_depth', 'MUNCHENB_sunshine', 'MUNCHENB_temp_mean', 'MUNCHENB_temp_min', 'MUNCHENB_temp_max', 'OSLO_cloud_cover', 'OSLO_wind_speed', 'OSLO_humidity', 'OSLO_pressure', 'OSLO_global_radiation', 'OSLO_precipitation', 'OSLO_snow_depth', 'OSLO_sunshine', 'OSLO_temp_mean', 'OSLO_temp_min', 'OSLO_temp_max', 'ROMA_cloud_cover', 'ROMA_wind_speed', 'ROMA_humidity', 'ROMA_pressure', 'ROMA_sunshine', 'ROMA_temp_mean', 'SONNBLICK_cloud_cover', 'SONNBLICK_wind_speed', 'SONNBLICK_humidity', 'SONNBLICK_pressure', 'SONNBLICK_global_radiation', 'SONNBLICK_precipitation', 'SONNBLICK_sunshine', 'SONNBLICK_temp_mean', 'SONNBLICK_temp_min', 'SONNBLICK_temp_max', 'STOCKHOLM_cloud_cover', 'STOCKHOLM_pressure', 'STOCKHOLM_global_radiation', 'STOCKHOLM_precipitation', 'STOCKHOLM_sunshine', 'STOCKHOLM_temp_mean', 'STOCKHOLM_temp_min', 'STOCKHOLM_temp_max', 'TOURS_wind_speed', 'TOURS_humidity', 'TOURS_pressure', 'TOURS_global_radiation', 'TOURS_precipitation', 'TOURS_temp_mean', 'TOURS_temp_min', 'TOURS_temp_max', 'VALENTIA_cloud_cover', 'VALENTIA_humidity', 'VALENTIA_pressure', 'VALENTIA_global_radiation', 'VALENTIA_precipitation', 'VALENTIA_snow_depth', 'VALENTIA_sunshine', 'VALENTIA_temp_mean', 'VALENTIA_temp_min', 'VALENTIA_temp_max', 'Date', 'Month']
In [4]:
# create a table that includes all columns for average temperature
temp = climate[[col for col in climate.columns if 'temp_mean' in col]]
temp.shape
Out[4]:
(22950, 18)
In [7]:
# view temp table column names
temp.columns
Out[7]:
Index(['BASEL_temp_mean', 'BELGRADE_temp_mean', 'BUDAPEST_temp_mean',
       'DEBILT_temp_mean', 'DUSSELDORF_temp_mean', 'GDANSK_temp_mean',
       'HEATHROW_temp_mean', 'KASSEL_temp_mean', 'LJUBLJANA_temp_mean',
       'MAASTRICHT_temp_mean', 'MADRID_temp_mean', 'MUNCHENB_temp_mean',
       'OSLO_temp_mean', 'ROMA_temp_mean', 'SONNBLICK_temp_mean',
       'STOCKHOLM_temp_mean', 'TOURS_temp_mean', 'VALENTIA_temp_mean'],
      dtype='object')
In [8]:
# create a box plot to view the temperatures
import seaborn as sns
sns.boxplot(data=temp)
plt.xticks(rotation=90)
Out[8]:
(array([ 0,  1,  2,  3,  4,  5,  6,  7,  8,  9, 10, 11, 12, 13, 14, 15, 16,
        17]),
 [Text(0, 0, 'BASEL_temp_mean'),
  Text(1, 0, 'BELGRADE_temp_mean'),
  Text(2, 0, 'BUDAPEST_temp_mean'),
  Text(3, 0, 'DEBILT_temp_mean'),
  Text(4, 0, 'DUSSELDORF_temp_mean'),
  Text(5, 0, 'GDANSK_temp_mean'),
  Text(6, 0, 'HEATHROW_temp_mean'),
  Text(7, 0, 'KASSEL_temp_mean'),
  Text(8, 0, 'LJUBLJANA_temp_mean'),
  Text(9, 0, 'MAASTRICHT_temp_mean'),
  Text(10, 0, 'MADRID_temp_mean'),
  Text(11, 0, 'MUNCHENB_temp_mean'),
  Text(12, 0, 'OSLO_temp_mean'),
  Text(13, 0, 'ROMA_temp_mean'),
  Text(14, 0, 'SONNBLICK_temp_mean'),
  Text(15, 0, 'STOCKHOLM_temp_mean'),
  Text(16, 0, 'TOURS_temp_mean'),
  Text(17, 0, 'VALENTIA_temp_mean')])
In [5]:
# add date column to temp table
temp_date = temp.copy()
temp_date['Date'] = climate['Date']
temp_date.head()
Out[5]:
BASEL_temp_mean BELGRADE_temp_mean BUDAPEST_temp_mean DEBILT_temp_mean DUSSELDORF_temp_mean GDANSK_temp_mean HEATHROW_temp_mean KASSEL_temp_mean LJUBLJANA_temp_mean MAASTRICHT_temp_mean MADRID_temp_mean MUNCHENB_temp_mean OSLO_temp_mean ROMA_temp_mean SONNBLICK_temp_mean STOCKHOLM_temp_mean TOURS_temp_mean VALENTIA_temp_mean Date
0 -0.528623 -1.016876 -1.099163 -0.114356 -0.105836 -0.927601 -0.106469 -0.182904 -1.370824 -0.097084 -0.988280 -0.265742 -0.186575 -1.280450 -0.124331 -0.391072 -0.257321 -0.668215 19600101
1 -0.582946 -1.107669 -1.110927 -0.367511 -0.370915 -0.825294 -0.892676 -0.212437 -1.043881 -0.232112 -0.691740 -0.353714 -0.368598 -0.539569 -0.650834 -0.415953 -0.335759 -0.548046 19600102
2 -0.257010 -1.084971 -1.063873 -0.509912 -0.532908 -0.940389 -0.490837 -0.389635 -0.741156 -0.487164 -0.853490 -0.403983 -0.550620 -0.876333 -0.650834 -0.615003 -0.210258 -0.067372 19600103
3 -0.555784 -1.209812 -1.146217 -0.525734 -0.577088 -1.042696 -0.316124 -0.493001 -0.910682 -0.472161 -0.624345 -0.642763 -0.417137 -0.775304 -0.943336 -0.764290 -0.069069 -0.998679 19600104
4 -1.003946 -1.209812 -1.087400 -0.320045 -0.444548 -0.978754 -0.403481 -0.552067 -0.862246 -0.307127 -0.381721 -0.906678 -0.332193 -0.926848 -0.621584 -0.503037 -0.037694 -1.509396 19600105
In [6]:
temp_date.to_csv(r'C:\Users\kacie\Box\Personal\CareerFoundry\Machine Learning\Climate Data\Avg Temp w Date.csv', index = False)

Madrid 1989 Analysis¶

Madrid 1989 Initial Analysis¶

In [10]:
# select only values for 1989 and drop the date column
temp_1989 = temp_date[temp_date['Date'].astype(str).str.contains('1989')]
temp_1989.drop(['Date'], axis = 1, inplace = True)
temp_1989.head()

C:\Users\kacie\AppData\Local\Temp\ipykernel_9656\224814569.py:3: SettingWithCopyWarning: 
A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
  temp_1989.drop(['Date'], axis = 1, inplace = True)
Out[10]:
BASEL_temp_mean BELGRADE_temp_mean BUDAPEST_temp_mean DEBILT_temp_mean DUSSELDORF_temp_mean GDANSK_temp_mean HEATHROW_temp_mean KASSEL_temp_mean LJUBLJANA_temp_mean MAASTRICHT_temp_mean MADRID_temp_mean MUNCHENB_temp_mean OSLO_temp_mean ROMA_temp_mean SONNBLICK_temp_mean STOCKHOLM_temp_mean TOURS_temp_mean VALENTIA_temp_mean
10593 -1.316301 -1.164415 -1.028583 -0.430800 -0.621268 -0.620680 -0.665550 -0.552067 -1.479805 -0.667200 -1.042197 -1.321402 -0.720508 -1.718244 -0.416832 -1.112626 -1.590770 -0.187540
10594 -1.452108 -1.334652 -1.169744 -0.683956 -0.724354 -0.927601 -0.717963 -0.625900 -1.407151 -0.667200 -1.230904 -0.956948 -0.805452 -1.600376 -0.577708 -0.615003 -1.700583 -0.277667
10595 -1.357043 -1.618381 -1.475592 -1.206090 -1.107247 -0.748564 -0.875205 -1.201794 -1.600895 -1.327334 -1.473528 -1.183161 -0.914666 -1.650891 -0.168206 -0.391072 -1.982960 -0.097414
10596 -1.465688 -1.607032 -1.628516 -1.095334 -0.812714 -1.145003 -0.560722 -1.157495 -1.818857 -0.952258 -1.635277 -1.333969 -0.841857 -1.886626 0.036546 -0.540359 -1.543707 -0.758341
10597 -1.058269 -1.958855 -1.828494 -0.826356 -0.842167 -0.889236 -0.875205 -0.876931 -1.721985 -0.862240 -1.473528 -1.019785 -0.671969 -1.836111 -0.797085 -0.453275 -0.978952 -0.548046
In [11]:
# create a 3D plot of weather station data for 1989 using plotly

import plotly.graph_objects as go

fig = go.Figure(data=[go.Surface(z=temp_1989.values)])
fig.update_layout(title='Weather Station Average Temperatures Over Days',
                  scene=dict(
        xaxis_title='Sation',
        yaxis_title='Day',
        zaxis_title='Average Temperature'
    ), 
                  autosize=True,
                  )
fig.show()
In [12]:
# create a df of days scaled
i = np.arange(0.01,3.66,0.01)
days = pd.DataFrame(data = i, columns = ['days'])
days
Out[12]:
days
0 0.01
1 0.02
2 0.03
3 0.04
4 0.05
... ...
360 3.61
361 3.62
362 3.63
363 3.64
364 3.65

365 rows × 1 columns

In [13]:
# find number of entries for 1989 table
n_rows = temp_1989.shape[0]
n_rows
Out[13]:
365
In [14]:
# create arrays for optimization function
X=days.to_numpy().reshape(n_rows,1)
#Represent x_0 as a vector of 1s for vector computation
ones = np.ones((n_rows,1))
X = np.concatenate((ones, X), axis=1)
y=temp_1989['MADRID_temp_mean'].to_numpy().reshape(n_rows,1)

X.shape, y.shape
Out[14]:
((365, 2), (365, 1))
In [15]:
# define the compute loss function
def compute_cost(X, y, theta=np.array([[0],[0]])):
    """Given covariate matrix X, the prediction results y and coefficients theta
    compute the loss"""
    
    m = len(y)
    J=0 # initialize loss to zero
    
    # reshape theta
    theta=theta.reshape(2,1)
    
    # calculate the hypothesis - y_hat
    h_x = np.dot(X,theta)
    #print(h_x)
    
    # subtract y from y_hat, square and sum
    error_term = sum((h_x - y)**2)
    
    # divide by twice the number of samples - standard practice.
    loss = error_term/(2*m)
    
    return loss
In [16]:
# define the gradient descent function
def gradient_descent(X, y, theta=np.array([[0],[0]]),
                    alpha=0.01, num_iterations=1500):
    """
    Solve for theta using Gradient Descent optimiztion technique. 
    Alpha is the learning rate
    """
    m = len(y)
    J_history = []
    theta0_history = []
    theta1_history = []
    theta = theta.reshape(2,1)
    
    for i in range(num_iterations):
        error = (np.dot(X, theta) - y)
        
        term0 = (alpha/m) * sum(error* X[:,0].reshape(m,1))
        term1 = (alpha/m) * sum(error* X[:,1].reshape(m,1))
        
        # update theta
        term_vector = np.array([[term0],[term1]])
        #print(term_vector)
        theta = theta - term_vector.reshape(2,1)
        
        # store history values
        theta0_history.append(theta[0].tolist()[0])
        theta1_history.append(theta[1].tolist()[0])
        J_history.append(compute_cost(X,y,theta).tolist()[0])
        
    return (theta, J_history, theta0_history, theta1_history)
In [18]:
%%time

# define starting values and run gradient descent function

num_iterations=10 
theta_init=np.array([[-10],[-5]]) 
alpha=0.1 
theta, J_history, theta0_history, theta1_history = gradient_descent(X,y, theta_init,
                                                                   alpha, num_iterations)

CPU times: total: 15.6 ms
Wall time: 12 ms
In [19]:
# show ending theta values
theta
Out[19]:
array([[-5.37579264],
       [ 2.41563456]])
In [20]:
# Show last 10 values for loss
J_history[-10:]
Out[20]:
[49.40629471118493,
 14.921604534545523,
 6.9884688413844245,
 5.058296927480842,
 4.489744585674478,
 4.23427832361378,
 4.054844643365118,
 3.8976801093281286,
 3.7504283805529344,
 3.6100942276962757]
In [21]:
# plot thetas and loss over iterations

fig, ax1 = plt.subplots()

# plot thetas over iterations
color='tab:blue'
ax1.plot(theta0_history, label='$\\theta_{0}$', linestyle='--', color=color)
ax1.plot(theta1_history, label='$\\theta_{1}$', linestyle='-', color=color)
# ax1.legend()
ax1.set_xlabel('Iterations'); ax1.set_ylabel('$\\theta$', color=color);
ax1.tick_params(axis='y', labelcolor=color)

# plot loss function over iterations
color='tab:red'
ax2 = ax1.twinx()
ax2.plot(J_history, label='Loss function', color=color)
ax2.set_title('Values of $\\theta$ and $J(\\theta)$ over iterations, Madrid 1989 Initial')
ax2.set_ylabel('Loss: $J(\\theta)$', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# ax2.legend();
fig.legend();
In [27]:
theta0_history
Out[27]:
[-8.07249075848817,
 -7.089733708473318,
 -6.55737838869459,
 -6.24045655936805,
 -6.027271132481126,
 -5.864709834727192,
 -5.727503232380554,
 -5.603614699199081,
 -5.48729703574407,
 -5.375792639018528]
In [28]:
theta1_history
Out[28]:
[-0.8907034571268202,
 1.0334824731171144,
 1.9198101150532585,
 2.3134943746541494,
 2.4736341912896505,
 2.523353129168908,
 2.521153177969456,
 2.4948253987799167,
 2.4575658380365755,
 2.41563455561203]
In [31]:
%%time

# define theta ranges and compute costs for range

# theta range
theta0_vals = np.linspace(-8.07,-5.2,100) #Look in the chart above for the limits of where theta0 and theta1 appear.
theta1_vals = np.linspace(-0.89,2.6,100) #Put those values as the first two "linspace" numbers in these lines
                                     
J_vals = np.zeros((len(theta0_vals), len(theta1_vals)))

# compute cost for each combination of theta
c1=0; c2=0
for i in theta0_vals:
    for j in theta1_vals:
        t = np.array([i, j])
        J_vals[c1][c2] = compute_cost(X, y, t.transpose()).tolist()[0]
        c2=c2+1
    c1=c1+1
    c2=0 # reinitialize to 0

CPU times: total: 2.47 s
Wall time: 2.5 s
In [32]:
# Plot Loss vs thetas

fig = go.Figure(data=[go.Surface(x=theta1_vals, y=theta0_vals, z=J_vals)])
fig.update_layout(title='Loss Function Across Thetas, Madrid 1989 Initial', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()
In [33]:
# Plot loss vs thetas with loss markers

line_marker = dict(color='#101010', width=2)
fig = go.Figure()
fig.add_surface(x=theta1_vals, y=theta0_vals, z=J_vals)
fig.add_scatter3d(x=theta1_history, y=theta0_history, z=J_history, line=line_marker, name='')

#The below line adds a graph of just the loss over iterations in a 2D plane
plt.plot(theta0_history, theta1_history, 'r+');
fig.update_layout(title='Loss Across Thetas, Madrid 1989 Initial', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()
In [34]:
#This plot shows the convergence similar to above, but only in the X/Y plane (there's no height)

plt.contour(theta0_vals, theta1_vals, J_vals, levels = np.logspace(0,10,1000))
plt.xlabel('$\\theta_{0}$'); plt.ylabel("$\\theta_{1}$")
plt.title("Contour Plot of Loss Across $\\theta$s, Madrid 1989");
plt.plot(theta0_history, theta1_history, 'r+');

Madrid 1989 Analysis Round 2¶

In [35]:
# rerun gradient descent function
num_iterations=200
theta_init=np.array([[-5],[2]])
alpha= 0.3
theta1, J_history1, theta0_history1, theta1_history1 = gradient_descent(X,y, theta_init,
                                                                   alpha, num_iterations)
In [36]:
# plot thetas and loss over iterations

fig, ax1 = plt.subplots()

# plot thetas over iterations
color='tab:blue'
ax1.plot(theta0_history1, label='$\\theta_{0}$', linestyle='--', color=color)
ax1.plot(theta1_history1, label='$\\theta_{1}$', linestyle='-', color=color)
# ax1.legend()
ax1.set_xlabel('Iterations'); ax1.set_ylabel('$\\theta$', color=color);
ax1.tick_params(axis='y', labelcolor=color)

# plot loss function over iterations
color='tab:red'
ax2 = ax1.twinx()
ax2.plot(J_history1, label='Loss function', color=color)
ax2.set_title('Values of $\\theta$ and $J(\\theta)$ Over Iterations, Madrid 1989 Rd 2')
ax2.set_ylabel('Loss: $J(\\theta)$', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# ax2.legend();
fig.legend();
In [37]:
theta1
Out[37]:
array([[-0.31757558],
       [ 0.24189401]])
In [38]:
# show last 10 values for loss history
J_history1[-10:]
Out[38]:
[0.4036592253315138,
 0.40365922532722304,
 0.40365922532345944,
 0.4036592253201582,
 0.40365922531726306,
 0.4036592253147233,
 0.4036592253124961,
 0.4036592253105426,
 0.40365922530882936,
 0.4036592253073262]
In [41]:
# define theta ranges and compute costs for range

# theta range
theta0_vals1 = np.linspace(-5.38,-0.318,100) #Look in the chart above for the limits of where theta0 and theta1 appear.
theta1_vals1 = np.linspace(2.42,0.242,100) #Put those values as the first two "linspace" numbers in these lines
                                      #Select with large margins, maybe +/- 10
J_vals1 = np.zeros((len(theta0_vals1), len(theta1_vals1)))

# compute cost for each combination of theta
c1=0; c2=0
for i in theta0_vals1:
    for j in theta1_vals1:
        t = np.array([i, j])
        J_vals1[c1][c2] = compute_cost(X, y, t.transpose()).tolist()[0]
        c2=c2+1
    c1=c1+1
    c2=0 # reinitialize to 0
In [42]:
# 3D plot of loss function vs thetas with loss markers

line_marker = dict(color='#101010', width=2)
fig = go.Figure()
fig.add_surface(x=theta1_vals1, y=theta0_vals1, z=J_vals1)
fig.add_scatter3d(x=theta1_history1, y=theta0_history1, z=J_history1, line=line_marker, name='')

#The below line adds a graph of just the loss over iterations in a 2D plane
plt.plot(theta0_history1, theta1_history1, 'r+');
fig.update_layout(title='Loss Function Across Thetas, Madrid 1989 Rd 2', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()

Oslo 1960 Analysis¶

Oslo 1960 Analysis Initial¶

In [43]:
# select only values for 1960 and drop the date column
temp_1960 = temp_date[temp_date['Date'].astype(str).str.contains('1960')]
temp_1960.drop(['Date'], axis = 1, inplace = True)
temp_1960.head()

C:\Users\kacie\AppData\Local\Temp\ipykernel_9656\357699861.py:3: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
Out[43]:
BASEL_temp_mean BELGRADE_temp_mean BUDAPEST_temp_mean DEBILT_temp_mean DUSSELDORF_temp_mean GDANSK_temp_mean HEATHROW_temp_mean KASSEL_temp_mean LJUBLJANA_temp_mean MAASTRICHT_temp_mean MADRID_temp_mean MUNCHENB_temp_mean OSLO_temp_mean ROMA_temp_mean SONNBLICK_temp_mean STOCKHOLM_temp_mean TOURS_temp_mean VALENTIA_temp_mean
0 -0.528623 -1.016876 -1.099163 -0.114356 -0.105836 -0.927601 -0.106469 -0.182904 -1.370824 -0.097084 -0.988280 -0.265742 -0.186575 -1.280450 -0.124331 -0.391072 -0.257321 -0.668215
1 -0.582946 -1.107669 -1.110927 -0.367511 -0.370915 -0.825294 -0.892676 -0.212437 -1.043881 -0.232112 -0.691740 -0.353714 -0.368598 -0.539569 -0.650834 -0.415953 -0.335759 -0.548046
2 -0.257010 -1.084971 -1.063873 -0.509912 -0.532908 -0.940389 -0.490837 -0.389635 -0.741156 -0.487164 -0.853490 -0.403983 -0.550620 -0.876333 -0.650834 -0.615003 -0.210258 -0.067372
3 -0.555784 -1.209812 -1.146217 -0.525734 -0.577088 -1.042696 -0.316124 -0.493001 -0.910682 -0.472161 -0.624345 -0.642763 -0.417137 -0.775304 -0.943336 -0.764290 -0.069069 -0.998679
4 -1.003946 -1.209812 -1.087400 -0.320045 -0.444548 -0.978754 -0.403481 -0.552067 -0.862246 -0.307127 -0.381721 -0.906678 -0.332193 -0.926848 -0.621584 -0.503037 -0.037694 -1.509396
In [47]:
# create a df of days scaled
i = np.arange(0.01,3.67,0.01)
days = pd.DataFrame(data = i, columns = ['days'])
days
Out[47]:
days
0 0.01
1 0.02
2 0.03
3 0.04
4 0.05
... ...
361 3.62
362 3.63
363 3.64
364 3.65
365 3.66

366 rows × 1 columns

In [48]:
# find number of entries for 1960 table
n_rows = temp_1960.shape[0]
n_rows
Out[48]:
366
In [49]:
# create arrays for optimization function
X=days.to_numpy().reshape(n_rows,1)
#Represent x_0 as a vector of 1s for vector computation
ones = np.ones((n_rows,1))
X = np.concatenate((ones, X), axis=1)
y=temp_1960['OSLO_temp_mean'].to_numpy().reshape(n_rows,1)

X.shape, y.shape
Out[49]:
((366, 2), (366, 1))
In [50]:
%%time

# define starting values and run gradient descent function
num_iterations=10 
theta_init=np.array([[-10],[-10]]) 
alpha=0.1 
theta2, J_history2, theta0_history2, theta1_history2 = gradient_descent(X,y, theta_init,
                                                                   alpha, num_iterations)

CPU times: total: 31.2 ms
Wall time: 14 ms
In [51]:
# show resulting thetas
theta2
Out[51]:
array([[-3.97546404],
       [ 1.70941059]])
In [52]:
# show last 10 values for loss history
J_history2[-10:]
Out[52]:
[104.56160896944326,
 25.419637216686407,
 7.656279054302163,
 3.617777963426695,
 2.6504504495448367,
 2.3724215732680713,
 2.2510818590203865,
 2.1672517132298346,
 2.0941832204597914,
 2.025797012726085]
In [53]:
# plot thetas and loss over iterations

fig, ax1 = plt.subplots()

# plot thetas over iterations
color='tab:blue'
ax1.plot(theta0_history2, label='$\\theta_{0}$', linestyle='--', color=color)
ax1.plot(theta1_history2, label='$\\theta_{1}$', linestyle='-', color=color)
# ax1.legend()
ax1.set_xlabel('Iterations'); ax1.set_ylabel('$\\theta$', color=color);
ax1.tick_params(axis='y', labelcolor=color)

# plot loss function over iterations
color='tab:red'
ax2 = ax1.twinx()
ax2.plot(J_history2, label='Loss function', color=color)
ax2.set_title('Values of $\\theta$ and $J(\\theta)$ Over Iterations, Oslo 1960 Initial')
ax2.set_ylabel('Loss: $J(\\theta)$', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# ax2.legend();
fig.legend();
In [54]:
theta0_history2
Out[54]:
[-7.172842210876891,
 -5.789402682355453,
 -5.089575426841161,
 -4.713943735247013,
 -4.492546543105828,
 -4.345002250855698,
 -4.233284936159303,
 -4.139395840802454,
 -4.054805922070791,
 -3.9754640371106436]
In [55]:
theta1_history2
Out[55]:
[-3.6730109989680537,
 -0.7015215158345369,
 0.6838345515699109,
 1.3196456158399232,
 1.6014613143518421,
 1.716298089985204,
 1.7525732282408233,
 1.7520842205541833,
 1.7345858112808796,
 1.709410592882017]
In [64]:
%%time

# define theta ranges and compute costs for range

# theta range
theta0_vals2 = np.linspace(-8,-3,100) #Look in the chart above for the limits of where theta0 and theta1 appear.
theta1_vals2 = np.linspace(-4,2,100) #Put those values as the first two "linspace" numbers in these lines
                                     
J_vals2 = np.zeros((len(theta0_vals2), len(theta1_vals2)))

# compute cost for each combination of theta
c1=0; c2=0
for i in theta0_vals2:
    for j in theta1_vals2:
        t = np.array([i, j])
        J_vals2[c1][c2] = compute_cost(X, y, t.transpose()).tolist()[0]
        c2=c2+1
    c1=c1+1
    c2=0 # reinitialize to 0

CPU times: total: 2.64 s
Wall time: 2.7 s
In [65]:
# 3D plot of loss function vs thetas with loss markers

line_marker = dict(color='#101010', width=2)
fig = go.Figure()
fig.add_surface(x=theta1_vals2, y=theta0_vals2, z=J_vals2)
fig.add_scatter3d(x=theta1_history2, y=theta0_history2, z=J_history2, line=line_marker, name='')

#The below line adds a graph of just the loss over iterations in a 2D plane
plt.plot(theta0_history2, theta1_history2, 'r+');
fig.update_layout(title='Loss Function Across Thetas, Oslo 1960 Initial', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()

Oslo 1960 Analysis Round 2¶

In [66]:
%%time

# define starting values and run gradient descent function
num_iterations=500
theta_init=np.array([[-3.97],[1.7]]) 
alpha=0.1 
theta3, J_history3, theta0_history3, theta1_history3 = gradient_descent(X,y, theta_init,
                                                                   alpha, num_iterations)

CPU times: total: 359 ms
Wall time: 361 ms
In [67]:
# plot thetas and loss over iterations

fig, ax1 = plt.subplots()

# plot thetas over iterations
color='tab:blue'
ax1.plot(theta0_history3, label='$\\theta_{0}$', linestyle='--', color=color)
ax1.plot(theta1_history3, label='$\\theta_{1}$', linestyle='-', color=color)
# ax1.legend()
ax1.set_xlabel('Iterations'); ax1.set_ylabel('$\\theta$', color=color);
ax1.tick_params(axis='y', labelcolor=color)

# plot loss function over iterations
color='tab:red'
ax2 = ax1.twinx()
ax2.plot(J_history3, label='Loss function', color=color)
ax2.set_title('Values of $\\theta$ and $J(\\theta)$ Over Iterations, Oslo 1960 Rd 2')
ax2.set_ylabel('Loss: $J(\\theta)$', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# ax2.legend();
fig.legend();
In [68]:
# show last 10 values for loss history
J_history3[-10:]
Out[68]:
[0.4732461066076224,
 0.4732461065593546,
 0.473246106513109,
 0.473246106468802,
 0.47324610642635095,
 0.4732461063856787,
 0.47324610634671094,
 0.47324610630937597,
 0.4732461062736051,
 0.47324610623933316]
In [69]:
# Show resulting thetas
theta3
Out[69]:
array([[-0.45432786],
       [ 0.20484415]])
In [70]:
%%time

# define theta ranges and compute costs for range

# theta range
theta0_vals3 = np.linspace(-3.97,-0.454,100) #Look in the chart above for the limits of where theta0 and theta1 appear.
theta1_vals3 = np.linspace(1.7,0.205,100) #Put those values as the first two "linspace" numbers in these lines
                                     
J_vals3 = np.zeros((len(theta0_vals3), len(theta1_vals3)))

# compute cost for each combination of theta
c1=0; c2=0
for i in theta0_vals3:
    for j in theta1_vals3:
        t = np.array([i, j])
        J_vals3[c1][c2] = compute_cost(X, y, t.transpose()).tolist()[0]
        c2=c2+1
    c1=c1+1
    c2=0 # reinitialize to 0

CPU times: total: 2.89 s
Wall time: 2.91 s
In [71]:
# 3D plot of loss function vs thetas with loss markers

line_marker = dict(color='#101010', width=2)
fig = go.Figure()
fig.add_surface(x=theta1_vals3, y=theta0_vals3, z=J_vals3)
fig.add_scatter3d(x=theta1_history3, y=theta0_history3, z=J_history3, line=line_marker, name='')

#The below line adds a graph of just the loss over iterations in a 2D plane
plt.plot(theta0_history3, theta1_history3, 'r+');
fig.update_layout(title='Loss Function Across Thetas, Oslo 1960 Rd 2', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()

Dusseldorf 2015 Analysis¶

Dusseldorf 2015 Analysis Initial¶

In [72]:
# select only values for 2015 and drop the date column
temp_2015 = temp_date[temp_date['Date'].astype(str).str.contains('2015')]
temp_2015.drop(['Date'], axis = 1, inplace = True)
temp_2015.head()

C:\Users\kacie\AppData\Local\Temp\ipykernel_9656\4248263576.py:3: SettingWithCopyWarning:


A value is trying to be set on a copy of a slice from a DataFrame

See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy
Out[72]:
BASEL_temp_mean BELGRADE_temp_mean BUDAPEST_temp_mean DEBILT_temp_mean DUSSELDORF_temp_mean GDANSK_temp_mean HEATHROW_temp_mean KASSEL_temp_mean LJUBLJANA_temp_mean MAASTRICHT_temp_mean MADRID_temp_mean MUNCHENB_temp_mean OSLO_temp_mean ROMA_temp_mean SONNBLICK_temp_mean STOCKHOLM_temp_mean TOURS_temp_mean VALENTIA_temp_mean
20089 -1.818785 -2.344726 -1.734387 -1.111156 -1.328146 -0.569527 -0.717963 -0.005705 -2.000492 -1.222313 -1.379174 -1.497345 -0.465676 -0.000746 -0.650834 -0.328869 -1.637832 -0.007287
20090 -1.615075 -1.629730 -1.357958 -0.430800 -0.753808 -0.416067 -0.455894 -0.005705 -1.370824 -0.757219 -1.271341 -1.019785 -0.174440 -0.000746 -0.036580 -0.204463 -1.198579 -0.007287
20091 -0.678010 -1.039574 -0.840368 -0.984579 -1.180880 -0.556739 -1.102331 -0.005705 -1.419260 -1.162301 -1.163509 -0.768437 -0.477811 -0.000746 -0.372957 -0.615003 -0.476948 -0.007287
20092 -0.922462 -1.141717 -1.052110 -0.968756 -1.180880 -0.633469 -1.451757 -0.005705 -0.922791 -1.162301 -0.974801 -0.944380 -0.878261 -0.000746 -1.031086 -0.888695 -0.649512 -0.007287
20093 -1.370624 -1.346001 -1.181507 -1.206090 -1.180880 -0.722987 -1.102331 -0.005705 -0.934900 -1.342338 -1.096113 -0.919246 -0.890396 -0.000746 -1.601465 -1.112626 -1.355455 -0.007287
In [73]:
# create a df of days scaled
i = np.arange(0.01,3.66,0.01)
days = pd.DataFrame(data = i, columns = ['days'])
days
Out[73]:
days
0 0.01
1 0.02
2 0.03
3 0.04
4 0.05
... ...
360 3.61
361 3.62
362 3.63
363 3.64
364 3.65

365 rows × 1 columns

In [74]:
# find number of entries for 2015 table
n_rows = temp_2015.shape[0]
n_rows
Out[74]:
365
In [75]:
# create arrays for optimization function
X=days.to_numpy().reshape(n_rows,1)
#Represent x_0 as a vector of 1s for vector computation
ones = np.ones((n_rows,1))
X = np.concatenate((ones, X), axis=1)
y=temp_2015['DUSSELDORF_temp_mean'].to_numpy().reshape(n_rows,1)

X.shape, y.shape
Out[75]:
((365, 2), (365, 1))
In [77]:
%%time

# define starting values and run gradient descent function
num_iterations=10 
theta_init=np.array([[-10],[-10]]) 
alpha=0.1 
theta4, J_history4, theta0_history4, theta1_history4 = gradient_descent(X,y, theta_init,
                                                                   alpha, num_iterations)

CPU times: total: 0 ns
Wall time: 8.98 ms
In [78]:
# plot thetas and loss over iterations

fig, ax1 = plt.subplots()

# plot thetas over iterations
color='tab:blue'
ax1.plot(theta0_history4, label='$\\theta_{0}$', linestyle='--', color=color)
ax1.plot(theta1_history4, label='$\\theta_{1}$', linestyle='-', color=color)
# ax1.legend()
ax1.set_xlabel('Iterations'); ax1.set_ylabel('$\\theta$', color=color);
ax1.tick_params(axis='y', labelcolor=color)

# plot loss function over iterations
color='tab:red'
ax2 = ax1.twinx()
ax2.plot(J_history4, label='Loss function', color=color)
ax2.set_title('Values of $\\theta$ and $J(\\theta)$ Over Iterations, Dusseldorf 2015 Initial')
ax2.set_ylabel('Loss: $J(\\theta)$', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# ax2.legend();
fig.legend();
In [79]:
theta0_history4
Out[79]:
[-7.160119639226236,
 -5.764821321174968,
 -5.057021614668866,
 -4.676898754606313,
 -4.453428666086943,
 -4.30531315859003,
 -4.193896377712347,
 -4.100788499969421,
 -4.017231461152261,
 -3.9390450464416045]
In [80]:
theta1_history4
Out[80]:
[-3.6579562478408993,
 -0.6635913312286617,
 0.7402167332136311,
 1.3885253973498446,
 1.6781842616696316,
 1.7977863156015375,
 1.8369514797927322,
 1.8382632347188155,
 1.8219513233788531,
 1.7976221183209533]
In [87]:
%%time

# define theta ranges and compute costs for range

# theta range
theta0_vals4 = np.linspace(-7.3,-3.7,100) #Look in the chart above for the limits of where theta0 and theta1 appear.
theta1_vals4 = np.linspace(-3.8,1.9,100) #Put those values as the first two "linspace" numbers in these lines
                                     
J_vals4 = np.zeros((len(theta0_vals4), len(theta1_vals4)))

# compute cost for each combination of theta
c1=0; c2=0
for i in theta0_vals4:
    for j in theta1_vals4:
        t = np.array([i, j])
        J_vals4[c1][c2] = compute_cost(X, y, t.transpose()).tolist()[0]
        c2=c2+1
    c1=c1+1
    
    c2=0 # reinitialize to 0

CPU times: total: 2.41 s
Wall time: 2.44 s
In [88]:
# 3D plot of loss function vs thetas with loss markers

line_marker = dict(color='#101010', width=2)
fig = go.Figure()
fig.add_surface(x=theta1_vals4, y=theta0_vals4, z=J_vals4)
fig.add_scatter3d(x=theta1_history4, y=theta0_history4, z=J_history4, line=line_marker, name='')

#The below line adds a graph of just the loss over iterations in a 2D plane
plt.plot(theta0_history4, theta1_history4, 'r+');
fig.update_layout(title='Loss Function Across Thetas, Dusseldorf 2015 Initial', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()
In [89]:
# show last 10 values for loss history
J_history4[-10:]
Out[89]:
[106.45434626797639,
 25.933426227859428,
 7.67937979748893,
 3.4917200309639043,
 2.483798436836207,
 2.1966502268911223,
 2.0747864056071386,
 1.9926403109897755,
 1.9217516194283522,
 1.855596563569945]

Dusseldorf 2016 Analysis Round 2¶

In [90]:
%%time

# define starting values and run gradient descent function
num_iterations=200 
theta_init=np.array([[-3.94],[1.80]]) 
alpha=0.2 
theta5, J_history5, theta0_history5, theta1_history5 = gradient_descent(X,y, theta_init,
                                                                   alpha, num_iterations)

CPU times: total: 156 ms
Wall time: 158 ms
In [91]:
# plot thetas and loss over iterations

fig, ax1 = plt.subplots()

# plot thetas over iterations
color='tab:blue'
ax1.plot(theta0_history5, label='$\\theta_{0}$', linestyle='--', color=color)
ax1.plot(theta1_history5, label='$\\theta_{1}$', linestyle='-', color=color)
# ax1.legend()
ax1.set_xlabel('Iterations'); ax1.set_ylabel('$\\theta$', color=color);
ax1.tick_params(axis='y', labelcolor=color)

# plot loss function over iterations
color='tab:red'
ax2 = ax1.twinx()
ax2.plot(J_history5, label='Loss function', color=color)
ax2.set_title('Values of $\\theta$ and $J(\\theta)$ Over Iterations, Dusseldorf Rd 2')
ax2.set_ylabel('Loss: $J(\\theta)$', color=color)
ax1.tick_params(axis='y', labelcolor=color)

# ax2.legend();
fig.legend();
In [92]:
# show last 10 values for loss history
J_history5[-10:]
Out[92]:
[0.3538808846209966,
 0.3538808762542618,
 0.35388086858059375,
 0.3538808615425821,
 0.3538808550875706,
 0.3538808491672675,
 0.353880843737378,
 0.35388083875727905,
 0.3538808341897115,
 0.35388083000050363]
In [93]:
# show resulting thetas
theta5
Out[93]:
array([[-0.47646428],
       [ 0.31428373]])
In [94]:
%%time

# define theta ranges and compute costs for range

# theta range
theta0_vals5 = np.linspace(-3.94,-0.476,100) #Look in the chart above for the limits of where theta0 and theta1 appear.
theta1_vals5 = np.linspace(1.80,0.314,100) #Put those values as the first two "linspace" numbers in these lines
                                     
J_vals5 = np.zeros((len(theta0_vals5), len(theta1_vals5)))

# compute cost for each combination of theta
c1=0; c2=0
for i in theta0_vals5:
    for j in theta1_vals5:
        t = np.array([i, j])
        J_vals5[c1][c2] = compute_cost(X, y, t.transpose()).tolist()[0]
        c2=c2+1
    c1=c1+1
    c2=0 # reinitialize to 0

CPU times: total: 2.41 s
Wall time: 2.42 s
In [95]:
# 3D plot of loss function vs thetas with loss markers

line_marker = dict(color='#101010', width=2)
fig = go.Figure()
fig.add_surface(x=theta1_vals5, y=theta0_vals5, z=J_vals5)
fig.add_scatter3d(x=theta1_history5, y=theta0_history5, z=J_history5, line=line_marker, name='')

#The below line adds a graph of just the loss over iterations in a 2D plane
plt.plot(theta0_history5, theta1_history5, 'r+');
fig.update_layout(title='Loss Function Across Thetas, Dusseldorf 2015 Rd 2', autosize=True,
                  width=600, height=600, xaxis_title='theta0', 
                 yaxis_title='theta1')
fig.show()

The optimization is not converging to 0 because the data is not linear but is being fit to a linear model. A polynomial should be used.¶

In [ ]: